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SUMMARY 

Design of supersonic combustors for scramjet engines requires signi- 
ficant detailed information. about the flow field. Traditional approaches 
based mainly on experimentation have been limited by difficulties in the 
acquisition of meaningful data. Recent advances in computational capa- 
bilities have made it possible to begin to predict complex flow fields in 
the combustor. Such analytical approaches could supplement the traditional 
approach by eliminating the need for some of the costly experimental para- 
meter studies. 

A three-dimensional parabolic flow program designed for supersonic 
combustors is evaluated to determine its capabilities. The mathematical 
foundation and numerical procedure are reviewed; simplifications are pointed 
out and commented upon. The program is then evaluated numerically by 
applying it to several subsonic and supersonic, turbulent, reacting and non- 
reacting flow problems. Computational results are compared with available 
experimental or other analytical data. Good agreements are obtained when 
the simplifications on which the program is based are justified. Limita- 
tions of the program and the needs for improvement and extension are pointed 
out. The three-dimensional parabolic flow program appears to be .potentially 
useful for the development of supersonic combustors. 
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National Research Council Senior Resident Research Associate, NASA Langley 
Research Center. 



INTRODUCTION 


Design of supersonic combustors for scramjet engines requires signi- 
ficant detailed information about the flow' field'. The traditional approach 
to combustor development has been.^mainly based on expensive, cut-and-try 
experimentation. Because of the highly turbulent and high-temperature 
environment in supersonic combustors, the acquisition of meaningful data 
is extremely difficult. Due to recent advances in computational techniques 
and computer capabilities, however, it is now possible to make computations 
for three-dimensional mixing-reacting turbulent flow fields in combustors. 
Analytical predictions for guiding the supersonic combustor development and 
experimentation will be very valuable to reduce costly hardware iterations 
based purely on experimentation. 

Several numerical approaches (for example, refs. 1-3) are available 
at present to predict three-dimensional turbulent mixing flow fields. 

Based on published results in the technical literature, the approach (ref. 1). 
developed by Spalding's group at Imperial College seems to be the most widely 
used. Satisfactory results have been obtained from their computer codes 
for many three-dimensional flow problems (for example, refs. 4-7). The 
three-dimensional parabolic flow computer program SHIP* (ref. 8) has been 
developed by Spalding's group for the NASA Langley Research Center and is 
used along with programs developed elsewhere for predicting combustor flow 
fields. This program is the subject of the present evaluation. 

SHIP (Supersonic Hydrogen Injection Program) is a later version of and 
supersedes HISS (Hydrogen Injection into a Supersonic Stream, ref. 9) 
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The general philosophy of the approach of Spalding's group has been 
reflected in many of their publications (for example, refs. 4,-10-13). 
That philosophy is to produce useful engineering results without being 
overly rigorous in formulations, and unduly time consuming or expensive. 
Therefore, in their approach, many assumptions and/or simplifications 
are used to achieve efficiency in both computer storage and computing 
time; however, some of them are difficult to justify rigorously. The 
accuracy and applicability of such an approach can only be established 
by testing the computations against experimental data or against compute 
tions produced by a more rigorously formulated approach. Some of the 
results of such an evaluation are given in this paper. 

The present paper is organized as follows. The fundamental formu- 
lation and numerical procedure are outlined and reviewed in this paper. 
Assumptions and/or simplifications are pointed out and commented upon. 
The computer program is then applied to several reacting and nonreacting 
turbulent and nonturbulent flow problems and results are compared with 
available experimental data or to other computations. Finally, the over 
all capabilities of the program are discussed. Additional development 
and needs for extensions related to the supersonic combustor flow field 
predictions are pointed out. 


SYMBOLS 

The International System of Units (SI) is used in this paper. 
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eonstants associated with tui?buTenee mode;! 
coefficient of shear stress, 
damping factor defined by eq., (5) 
injector diameter . 
mass concentration 

total enthalpy or shape factor of turbulent boundary layer, 

6*/e 

turbulence kinetic energy 

turbulence length scale 

mixing length 

Mach number 

pressure 

Prandtl number 

Reynolds number 

source term, eq. (1) 

jet spacing 

Schmidt number 

temperature 

velocity components 

free-stream velocities of two mixing streams, eq. (6) 

nondimensional velocity of turbulent boundary layer, eq. (3,) 
rectangular coordinates 

nondimefisipnal cpordinate in turbulent boundary layer, eq. {3} 
Injection location 
exchange coefficient 
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Subscripts: 

eff 

E, W, S, N 

Z 

0 

P 
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boundary layer or mixing layer thickness 
displacement thickness 
energy thickness 

turbulence dissipation energy rate 

momentum thickness 

constant, eq. (4) 

constant, eq. (4) 

viscosity 

nondimensional coordinate in mixing layer, eq. (6) 
densi ty 

spreading constant of mixing layer 
shear stress 

general dependent variable, eq. (1) 
effective 

east, west, south, and north, respectively 

laminar 

total 

particular node or pitot condition 
turbulent 


w 


wall 

general dependent variable 



REVIEW OF THE COMPUTATIONAL APPROACH 


In this section the three-dimensional parabolic flow computer program, 
SHIP, developed for the NASA Langley Research Center (ref. 8), is intro- 
duced. Its theoretical foundation, numerical procedure, and simplifica- 
tions are discussed. 

The computer program is designed for calculating three-dimensional, 
turbulent, reacting, parabolic flows (either external or internal). The 
flow field considered in the program is a flow confined in a rectangular 
parallelepiped; any of the four lateral boundaries can be a wall, a symmetry 
plane, or a free-stream condition along a given surface. For walls, the 
distance of each wall from a reference plane may be specified as an arbi- 
trary smooth function of distance along the main flow direction. The 
main flow can be either subsonic or supersonic. Variable specific heats 
are used for different species in the flow and the mixture satisfies the 
equation of state of a perfect gas; four equilibrium chemical reactions are 
al 1 owed . 

The SHIP program was developed based on the Eulerian formulation in 
a rectangular coordinate system (x, y, z) with the z-axis in the main flow 
direction. The mean flow velocity components (u, v, w) , pressure (p), 
density (p), total enthalpy (H), and hydrogen mass concentration (f) of 
the three-dimensional turbulent mixing-reacting flow are governed approxi- 
'mately by the Navier-Stokes equations together with the equation of state 
for perfect gases and a- specie equation. To account for the turbulence, 
the laminar (molecular) viscosity (y^) is replaced by an effective viscosity 
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laminar (molecular) Prandtl number Pr^ and Schmidt 

number Sc^ are replaced by their empirical effective values to be discussed 

later. The turbulent viscosity (y^) is determined by way of a "k-e" two- 

equation turbulence model. Analogous to the molecular viscosity defined 

in the kinetic theory, the turbulent viscosity is related to the density, 

a turbulence velocity scale (e.g. square-root of the turbulence kinetic 

energy k) and a turbulence length scale l. At high Reynolds numbers, z 

3/2 

is proportional to k /e, with e the turbulence dissipation energy rate. 

2 

Hence p k /e with being an empirical constant, and k and e are 

determined by a set of transport equations with several additional empirical 
constants (ref. 14). 

One of the important simplifications for the development of the pro- 
grams is the parabolic flow assumption (ref. 4). When there exists a 
predominant flow direction, when the diffusion of mass, momentum, energy, 
etc., can be neglected in that direction, and when the dov^^nstrealn pres- 
sure field has little effect on the upstream flow field, the coordinate 
in the main flow direction (z) becomes a "one-way" coordinate. The set 


of governing equations reduces mathematically to the parabolic type, and 
can be numerically solved in succeeding cross-stream (x-y) planes in 
the main-stream direction. Because of this simplification, a three-dimen- 
sional program requires only two-dimensional computer storage. Consequently, 
computer storage and running time are greatly reduced. It should be' 
noted that the classical Prandtl boundary layer flow is a two-dimensional 
parabolic flow, but, unlike the simplification used here, the Prandtl 
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boundary layer equations can be justified rigorously by an order of magni- 
tude or a perturbation estimation. 

The simplified governing equations for the present parabolic flow can 
be represented by the following general differential form, 

■|t (‘"''t' ■ ^4. if' T 7 ■ '■4. if' ■ S 


where 9 is a general dependent variable and is a general exchange 
coefficient. When 0 = 1, u, v,w, H, f, k, ore, equation (1) corresponds, 
respectively, to the continuity, three components of momentum, energy, 
concentration, turbulence kinetic energy, or turbulence dissipation energy 
rate equation. On the left-hand side of equation (1), the first term 
represents convection in the main flow (z) direction, and the second and 
third terms represent, respectively , the sum of convection and diffusion 
in the lateral x and y directions. The term S, on the right-hand side 
is a source term which includes all other terms left from each of the 
simplified differential equations, such as the pressure gradient terms 
in the three momentum equations. The appropriate exchange coefficients 
and source terms for variables 4 ) corresponding to different 
conservation equations are listed in Table I.- (Note that some of the 
source' terms are different from those' in refs. 8 and 9, and Have been im- 
proved during the present evaluation.) The general effective exchang'd 
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• coefficient is composed of two parts, a turbulent and a laminar; i.e. 


^eff 

‘"'"eff,(|) 






where Pr. , and Pr. , are, respectively, the turbulent and laminar 
t,9 X,(p 

Prandtl numbers. The values of Pr^, Pr^ and the constants, C-j and 
in Table I are usually determined empirically (for example, ref. 14). 

The numerical formulation is based on a finite-difference form of 
equation (1). A "staggered" grid system is used in the x-y plane (ref. 8). 
The advantage of such a system is that velocity components u and v are 
stored just at the point (midv/ay between the two neighboring nodes) at 
which they are needed for the calculation of the convection, and the 
pressures are stored so as to make it easier to calculate the pressure 
gradients for calculating u and v. Thus, control volumes at each node 
are different for u and v as compared with w and the other (j>'s. By taking 
volume integrals of equation (1) over respective control volumes, a set 
of difference equations can be obtained. The volume integrations of the 
terms 5— (pwcj)) and S, in equation (1) are performed by assuming that 

o 2 . $ 

the values of (j), p, and w at a node point P are constant over the entire 

g 

control volume. The term -r — contains the difference of values at stations 

o 2 

z and z + Az. The volume integrals of the other two terms in equation 
(1) give rise to the surface integrals of the convective and diffusive 
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fluxes across the boundaries of tho control volume., A' proper representa- 
tion of these terms is essential to the convergence of numerical compu- 
tation. To provide numerical convergence and accuracy, a "hybrid" scheme is, 
used. This scheme is a combination of central and upwind differences; 
the rationale of this scheme and some experience with its use are described 
in references 10 and 15. When the integrations of various terms in equa- 
tion (1) are expressed in the manner described above, the following general 
form of the difference equations is obtained at an arbitrary node P, 


♦-P ' ''n‘‘n 




''E'f'E 


'+ Vw 


+ B 


( 2 ) 


Here, the values of 4>'s pertain to station z + Az and the values of A's 
and B pertain to station z. The subscripts N, S, E, and W denote, respec- 
tively, the neighboring north, south, east, and west nodes. The deriva- 
tion of equation (2) from equation (1) is given in the Appendix. 

In the present program, the set of difference equations (2) together 
with other auxiliary equations are not solved point-by-point simultane- 
ously, but by a so-called SIMPLE (for Semi-Implicit Method for Pressure 
Linked Equations) method (ref. 4). Although the flow problem is highly 
nonlinear, an economical noniterative procedure is' followed. The three 
■ velocity components are -solved -from -their respective (<(> = u,- v, w) difference 
equations in terms of a guessed pressure field. This guessed pressure field 
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1s then corrected by the continuity equation, i.e. equation (1) with ({) = 1. 
The pressure correction satisfies a difference equation similar to equation 
(2). After the guessed pressure field and three components of velocity 
field have been corrected by the pressure correction, the difference 
equations (2) for H,.f, k and e are solved sequentially. Temperature, 
density, mass concentrations of species and other auxiliary quantities 
are determined noniteratively by their appropriate relations. 

By knowing appropriate boundary conditions, each difference equation 
of equation (2) can be solved readily by several computer algorithms. 
However, it is usually adopted by the Spalding group that equation (2) 
is solved in a line-by-line iterative fashion by the successive use of a 
standard tri -diagonal matrix algorithm in the x and y directions. Great 
economy of computer time has been found by using this procedure (ref. 11). 

One additional feature of the programs is the handling of boundary 
conditions for these difference equations. The boundary conditions are 
specified by the values of appropriate fluxes across the boundaries. For 
a free or symmetry boundary, fluxes are automatically set to zero. For 
solid walls, to avoid using a large number of grid points to compute steep 
gradients in the sublayers, simplifications are employed. On such bound- 
aries, the v/ell- known (turbulent boundary layer) wall functions are used 
at near-boundary grid points located in the turbulent flow region. 

Moreover, the wall shear in the wall functions is replaced by the local 
Reynolds stress which, in turn, is related' to the local turbulence kinetic 
energy in an equilibrium high Reynolds number turbulence flow. Near 
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the. intersections of two walls, the corner fliows.. are not specially 
treated.. 

In the SHIP Program, the wall boundaries are aTTowed to vary along 
the main flow direction; the coordinates chosen are no longer orthogonal. 
However, it is stipulated that two of the coordinates (x, y) are maintained 
mutually orthogonal throughout the flow field, while the third (z) axis 
is permitted to depart from orthogonality with' respect to the other two. 

The error thus introduced is of the order of the tangent of the angle of 
departure from orthogonality. 

As reviewed above, the capabilities and limitations of the present 
parabolic flow programs depend primarily on the simplifications introduced 
and on the turbulence modeling. At present, it is believed that the 
"k-e" two-equation turbulence model with appropriate constants is suitable 
for high Reynolds number flows. In low Reynolds number flow regions, 
such as near a v/all or in a transition region, the turbulence model 
requires some special attention. 

Simplifications related to the development of the present program 
can be generally summarized and classified as follows: 

Simplifications which limit the range of the application of the 
program . - They are: (1) the parabolic flow assumption; and (2) the non- 
orthogonality of coordinates. The former excludes computations of flows 
with recirculation and/or with large pressure effects from downstream. 

The latter excludes computations with large variations of wall geometries. 
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Simplifications which affect the accuracy of the .program . - They 
are the noniterative SIMPLE procedures and the iterative procedure of 
solving difference equations. The accuracy of the program due to these 
simplifications may be improved by increasing grid numbers, using smaller 
forward steps and/or increasing numbers of iterations. 

Simplifications which mav affect the local flow field . - One 
such simplification is the application of wall -functions at near-boundary 
grid points for solid boundaries. As described previously, the wall shear 
in the wall functions is approximated by using the local Reynolds stress 
which, in turn, is related to the local turbulence kinetic energy in an 
equilibrium high Reynolds number turbulent flow. For practical applica- 
tions, the local Reynolds stresses at the near boundary grid points are not 
equal to the wall shear if the points are too far from the wall boundary; 
on the other hand, the high Reynolds number relationship between the local 
Reynolds stress and turbulence kinetic energy is not appropriate at the 
near boundary grid points If they are too close to the wall {ref. 14). 

Hence, using such a simplification in the v/all function, the program may not 
yield accurate wall shear stresses and heat transfer. Similarly, using 
such a simplified wall boundary condition, the program may not predict 
accurate flow fields near the wall. However, this flow field may be im- 
proved by using correct wall shear and/or by applying low Reynolds number 
corrections in the wall functions. 

The last simplification is hot related to the development of the 
programs, but is in the application. That is how to specify the conditions 
at the initial station. Mean flow quantities, such as velocity components, 



pressure and density, are usually measured and relatively easy to specify. 
The initial i;urfaulence kinetic energy and dissipation energy rate are • 
usually not measured and must be estimated for each computation. Currently, 
the turbulence kinetic energy and dissipation energy rate are estimated 
by means of the mixing-length hypothesis. The effects of such an estimate 
on some flow fields are significant and will be given in the next section 
along with comparisons of computational results and experimental data. 

APPLICATION OF THE COMPUTER PROGRAM 

The parabolic flow computer program SHIP has been applied to several 
flow fields. The objectives of such applications are: (!) to test the ' 

accuracy against available experimental data or other computational re- 
sults; (2) to establish or verify the values of empirical constants associ- 
ated with the "k-e“ two-equation turbulence model; and (3) to determine 
the limitations of this .program. Numerical computations were performed 
for the turbulent boundary layer flow of Wieghardt (ref. 16), the tur- 
bulent mixing layer flow of Brown and Roshko (ref. 18), compression and 
expansion of supersonic flow in a two-dimensional duct,. mixing of a jet 
normal to a supersonic stream (ref. 22), and supersonic combustion and 
mixing of a hydrogen wall jet in a duct (refs. 23, 24). 

Turbulent Boundary Layer Flow 

The two-dimensional turbulent boundary layer over a flat plate is 
perhaps the simplest and most-examined turbulent flow in the presence 
of a solid boundary. Such a flow satisfies the parabolic flow assumption. 
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and the program should be applicable. Wieghardt's flat plate flov.' examined 
in the 1968 Stanford Conference {ref. 16) is chosen as typical of the stan- 
dard flat-plate boundary layer. The free-stream fl6\^ velocity was 33 m/sec 
and the model was a waxed-plywood plate with a blunt leading edge fitted 
with a small trip wire. The mean-velocity profiles were measured by a probe 
rake at several downstream stations. The first profile at z = 0.087 m is 
probably at about the minimum Reynolds number for turbulent flow. The 

tunnel turbulence level was about 0.25 percent and an average kinematic 

2 

viscosity was about 0.151 cm /sec. 

Since the present program cannot predict boundary layer transition 
and the present two-equation turbulence model is not appropriate in low 
Reynolds number flows, computations by the program are started at 
z = 0.187 m (Re ^ 4.1 x 10^) and Wieghardt's measurements at that station 
are used for the initial conditions. Pressure, density, and temperature 
are assumed to be constant across the boundary layer and equal to their 
free-stream values. The ten initial velocity measurements across the boun- 
dary layer are matched by the following wall function near the wall 
(y'*’ < 50): 

+ + 
w = y 

w^ = 5.5 + 2.5 £ny^ 

Here w^ = — and 

yf\/p~ 


4 “ 


for 0<y <11.4 


for 11 .4 < y <50 


(3a) 

(3b) 


_ yyf- 


T P 




are, respectively, the nondimensional 
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velocity and coordinate perpendicular to the walU is :the wall shear. 

The initial turbulence kinetic energy and dissipation energy rate are ' 
estimated from the initial velocity profile by the- mixing length hypbthesis 
(ref. 17). The mixing-length is defined as follows, 

= X6 for y > I 6 (4a) 

where 

D = 1 ~ exp [-y(xp)^'^^/y^A] (5) 

is a damping factor. 6 = 0.005m is the boundary thickness at the station 
z = 0.187 m and the empirical constants X = 0.09, k = 0.435, and A = 26* 
Computations were performed to downstream stations with empirical 
Spalding constants associated with the "k-e'‘- two equation turbulence model 
(ref. 8): = 1.44, = 1.92, Cp = 0.09, laminar Prahdtl number Pr^^ = 0.7 

and turbulent Prandtl numbers Pr^ ^ - 1.0 for <J> = v, w, M and k, arid 1.3 for <j) 
At the near boundary grid points located jUst outside the sublayers a 
constant wall shear coefficient (at z = 0.187 m) = 0.00424 was used 
in the wall function. Results of velocity profiles at z = 0.187,' 0.287, 
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0.387, 0.487, and 0.637 m are compared with Wieghardt's data-in figure 1. 

The integral parameters, displacement thickness 6*, momentum thickness 0, 
energy thickness 6** and the shape factor H (= 5^6) compared with 
Wieghardt's data in figure 2. Comparisons of computations with experiments 
in figures 1 and 2 are in excellent agreement. 

Computations were also performed with the same set of empirical con- 
stants but with Reynolds stresses replacing the wall shear in the wall 
function. A typical velocity comparison (at z = 0.487 m) is shown in 
figure 3. The agreement with experimental data is not as good as that of 
the previous calculation. Due to the disagreement in velocity profiles, 
the comparisons of the integral parameters are also bad. This kind of 
disagreement should be expected, since the development of the turbulence 
is due to the presence of the wall; i.e. the v/all boundary conditions 
dominate the boundary layer flow field. This comparison sustains the 
conment made about the wall boundary condition in the previous section. 

Turbulent Mixing Layer Flow 

The two-dimensional turbulent mixing layer is the simplest turbulent 
mixing flow without influence from a wall. Recent experiments by 
Brown and Roshko (ref. 18) are used as the present test cases. The objec- 
tive is to compare the mixing profiles of the present computation with 
those of the experiments. 

Brown and Roshko 's experiments were conducted for mixing of two streams * 
separated upstream by a splitter plate. The streams having different 
velocities but the same density or different densities were mixed. Profiles 
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of velocities and densities were reported at several downstream' stations 
from the trailing edge of the splitter plate. Since- a trailing edge flow 
is theoretically not a parabolic flow, the present computations v/ere then 
started from the first station of their measurements. The initial turbu- 
lence kinetic energy and dissipation energy rate were estimated by the 
mixing-length hypothesis (ref. 17), The mixing-length is related 
to the thickness of the mixing layer 6, which, at a station, z, is the 
distance between two streamlines y-^ and y^ defined by [W.j-w(y.| )]/{W.j-W 2 } = G.l. 
and [W.j-w(y 2 .)]/(W.j-W 2 ) = 0.9, respectively. For twQ-dimens.io.nal plane flows, 

= 0.125 6 (ref. 19). 

Mixing of the two airstreams with different velocities (10 and 1.43 
m/sec) is considered first. The measured velocity data at the first 
station (z = 0,0254 m) and the initial profile for the computation are 
shown in figure 4; the thickness of the mixing layer at this station is 
0.0038 m. Since the distance of this station from splitter plate is about 
1000 tines the momentum thickness of the boundary layer leaving the plate, 
the velocity profile has already exhibited the similar nature of the 

free mixing layer. As noted in ref. 18, the virtual origin (z^) of the 

-3 

similar mixing layer is about 6.35 x 10 m upstream of the trailing edge. 

In the computation, the values of following empirical constants associated 
with the *‘k-e" two-equation turbulence model were used (ref. 19): C.j = 1.44, 

C 2 = 1.92, Cp = 0.09, Pr^ = 0.7, Pr^ ^ = 1.0 for ^ = v, w, and k, 1.3 for (f> = e and 
0;'5 for (})■= H, Pressure, density, and temperature were taken to be constants 
across the mixing layer at the initial station. 
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Computational results of velocity at downstream stations (2 = 0.0508, 
0.0635, 0-0762 and 0.0889 m) are compared with the experiments in figure 5. 
Comparisons are in good agreement except some discrepancies do occur at the 
edges of mixing layer where the experimental data deviate from the compu- 
tational results possibly due to large structures of the turbulence. A 
computation was also performed by fitting the initial velocity with a 
similar profile (ref. 20), 


W, + W, - VL 

w= -^^[1 + w^erf C] 


( 6 ) 


where W-j = 10 m/sec, = 1.43 m/sec; erf (C) is the error function with 
the argument 5 and ct = 13.5, y^ = -2.9 x 10~^m. The re- 

sults are indistinguishable from the previous computations, so the compu- 
tations are not shown separately in figure 5. 

The second test case is the mixing of flows with different media (N2 
and He) and with different velocities (10 and 3.78 m/sec, respectively, 
for N2 and He). The velocity and density data at the first station (z = 
0.0254 m) and the initial profiles for the computation are shown in figure 6 
(with = -0.0254 m). The thickness of the mixing layer at this station 
was 0.0031 m. In this computation, the same empirical constants as the air- 
air mixing case were used; in addition, the turbulent Prandtl number Pr^ = 
0.5 for the concentration was used in the specie equation. Computational 
results are compared with the corresponding experimental data at several 
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downstream' stations (z = 0.0508, 0.0762 and 0.10,16 m)- in figure 7. Agree- 
ments are good for the- velocity profiles, but not as good for the density 
profiles. 

The* discrepancies 'in the -density comparison led to. the calculation of- 
another mixing layer problem; i.e;, and He with velocities 1.43 and 10 
m/sec, respectively. Again, the velocity and density data at z = 0.0254 m 
were used as the initial conditions; they are shown in figure 8.. Since the 
thickness of the velocity and density mixing layers are not equal, the al- 
gebraic mean thickness was used (6 = 0.0064 m). Computational results of 
velocity and density are plotted in figure 9 to compare with experimental 
data at different downstream stations (z = 0.0508, 0.0635, 0.0762 and 
0.0889 m). Again, the agreements are good for the velocities; the com- 
parisons of densities are not good, especially in the part of mixing layer 
with higher density. 

Computations of these mixing layers have shown that the estimate of 
initial turbulence kinetic energy and dissipation energy rate has a signi- 
ficant effect on the downstream mixing profiles. For mixing of streams 
with equal densities, the Prandt.l mixing- length hypothesis with a constant 
mixing length seems to be appropriate for such an estimate. For mixing of 
streams with different densities, however, such an estimate together with 
the "k-e" two-equation turbulence model does not predict the same similar 
nature of the mixing, layer as. the experiment did and, consequently, yields . 
discrepancies. Different comparisons can- be obtained by trying different 
estimates of the initial turbulence kinetic energy. For example, by 
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assuming a mixing-length which varies with density (i.e. , S. = p^6) 

fn 97b 

in the last test case, a better agreement with the experiment was 
obtained. A typical comparison of such results at z = 0.0889 m is shown 
in figure 10. 

At present, there are no general methods which are also simple (like 
the Prandtl mixing-length hypothesis) to estimate "correct" initial turbu- 
lence kinetic energy and dissipation energy rate for use with the "k-e" 
two-equation turbulence model in a medium with large density variations. 

For the present mixing layer problem, two approaches may be possible to 
provide the initial data of the turbulence kinetic energy and dissipation 
energy rate. First, since the mixing layers at and downstream of the first 
station of measurements are already similar, a similarity consideration could 
be used to reduce the set of governing (partial differential) equations to 
ordinary differential equations. A "correct" and also similar, turbulence 
kinetic energy or dissipation energy rate profile could then be related to 
other known flow variable profiles and their derivatives at that station. 

Second, the "correct" turbulence kinetic energy (and other flow variables) 
as the initial conditions to the present parabolic flow program may be 
obtained by directly solving an elliptic- type turbulent-trai ling-edge flow 
problem. Since the mixing-length hypothesis is appropriate for the medium 
on either side of the splitter plate, the elliptic solution should be free 
of uncertainties to determine the upstream boundary conditions of the turbulence 
kinetic energy and dissipation energy rate. - 
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Compression and Expansion of Supersonic Flow 

In this section, supersonic flow in a two-dimensional channel with 
variable heights is computed by the SHIP Program. The objective is to 
demonstrate how compressions and expansions are handled by the parabolic 
flow computer program. 

The geometry of a two-dimensional channel, and the wave pattern in 
the channel are sketched in figure 11. The upper and lower walls of the 
channel are parallel and 0.2 m apart initially. At z = 0.1 m, the lower 
v/all is deflected by an angle a(= 3.366°) into the flow; at z = 0.3 m, the 
lower wall is deflected away from the flow by the same angle a. The flow 
is uniform initially with M = 1.5, p = 0.1 MN/m^, and T = 294 K. Computa- 
tions are performed for two cases: one with zero laminar viscosity and 

one with a variable laminar viscosity. In both cases, however, the turbu- 

-30 

lence viscosity and kinetic energy are set small (y.j. ~ 10 ) and 

kept constant. No-slip boundary conditions at walls are used in 
both cases. 

Corputational results of pressures at the upper and lower walls are 
plotted in figure 12. For comparison, the pressures computed by a tv/o- 
dimensional inviscid shock fitting method (ref. 21) are also plotted*. 

It may be noted that the shock fitting method employs a finite-difference 
numerical scheme; however, the locations and strengths of the shock dis- 
continuities are treated explicitly. By comparing the computational results 

*Acknowl edgment is made to J. P. Drummond foT his help in providing data by 
use of Program SEAGULL. 
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of two' methods, the results of the SHIP program show delays in pressure 
responses at both shocks and expansion corners. This kind of discrepancy 
is due to the finite-difference numerical scheme and the parabolic flow 
assumption. The former causes the smear of shocks and expansions at corners, 
whereas the latter prohibits transmission of disturbances to upstream 
points even in the subsonic flow region near the wall. Since a distur- 
bance can transmit laterally (and to downstream) by diffusion and con- 
vection, pressures in the laminar case show slightly better comparisons 
with those of the shock fitting method than those of the zero-viscosity 
case in the compression regions. 

Pressure variations across the channel at three stations (z = 0.24, 

0.35 and 0.48 m) are shown in figure 13. Comparisons of the two present 
computations with the shock fitting method are reasonably good except 
directly behind a shock or at an expansion corner due to the reasons stated 
above. 

Comparisons in figures 12 and 13 also show that a better agreement is 
obtained in regions with gradual pressure variations than with sharp changes. 
In the actual turbulent flow near a wall, shocks are smeared by large 
diffusive effects, so the pressure changes are usually more gradual than in 
the inviscid flow. Therefore, it is believed that the prediction of the 
present program will be in better agreement with experiments in a supersonic 
turbulent channel flow than with the results of an inviscid flow- 
field calculation. 

The computational results discussed above were obtained by using 60 
grid points equally spaced across the two-dimensional channel and the 
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forward step size being equal to 1/40 of the channel height. By varying 
the number of grid points and the step size, it was found that the numeri- 
cal results do not have appreciable changes for grid points over 30 and for 
the step size less than 1/40 of the channel height. However, in the in- 
viscid computations, the steep pressure variation across shocks can be 
improved by taking smaller forward step size. 

Mixing of Jet Normal to a Supersonic Stream 

In this section, the parabolic flow computer program is applied to 
a case of three-dimensional mixing of cold hydrogen injected normal to a 
supersonic airstream. Such experiments were reported in reference 22. 

The flow field and the arrangement of injectors at the injection station 
(z.) are sketched in figure 14. Experimental surveys of flow field were 

\J 

made at several downstream stations ({z-z.)/d = 7, 30, 60, 120, and 240 

J 

-3 

with the injector diameter d = 1.026 x 10 m) for two injector-arrange- 
ments (s/d = 6.25 and 12.5). 

Recognizing the occurrence of recirculating flow regions directly 
downstream of the injectors, the parabolic flow program was used to compute 
the flow field far downstream of the injectors where the recirculation 
is not present. The y-z plane was centered on the central jet; the two 
mid-planes (^ =+ 3.125 for the case of s/d = 6.25) between the central 
jet arid its neighbors, were chosen as the two "symnetric" boundaries. The 
other two boundaries were the wall (y = 0) and a free-stream boundary 
(y = H = 0.1 m). The experimental data (ref. 22) at the station 30 in- 
jector diameters downstream of the injection were used as initial conditions. 
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Since measurements were taken mostly along the symmetric plane (x = 0} 
at that station, very few 'measurements are available, at other locations 
at that station. In order' to make the experimental data fit onto the 
initial data (x,y) plane with 12 x 25 grid points, some flow quantities 
were assumed at several additional locations including the boundaries. 

By using a cubic spline interpolation routine, initial data at the 12 x 25 
grid points were obtained. Typical initial profiles of w, T, f and p 
at x/d = 0 and 2.19 are shovm in figure 15; the corresponding experimental 
data at x/d = 0 are plotted for comparison. 

Because of the uncertainty due to the large number of assumed initial 
values, the computation was not intended for detailed comparisons of the 
accuracy with experiments. Thus, for simplicity, the initial turbulence 
kinetic energy and dissipation energy rate v/ere assumed to be zero, and 
the initial lateral velocity components were also taken to be zero. Com- 
putations were then carried to the downstream stations. Typical profiles 
of w, T and f at x/d = 0 and 2.19 are shown in figures 16 and 17 for the 
stations (z-z.)/d = 60 and 120, respectively. The available experimental 

V 

data at x/d = 0 are also plotted. The comparisons show a general qualita- 
tive agreement between the experiments and computations. 

A similar computation for the case of s/d = 12.5 was also performed. 
Typical results of w, T and f at x/d = 0 and 5,0 and their comparisons 
with corresponding experiments at x/d = 0 are shown in figure 18. Again, 
the computations are in qualitative agreement with the experiments. To 
make detailed comparisons of three-dimensional computations with experiments. 
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a complete survey of flow fields is necessary and sliould be carefully 
planned. 


Supersonic Combustion and Mixing of Hydrogen in a Duct 

Most of turbulent reacting (combustion) experiments were conducted 
for axi symmetric configurations which are not appropriate for the direct 
application of the present program because of the rectangular coordinate 
system. Burrows and Kurkov (refs. 23, 24) performed supersonic combustion 
of hydrogen and turbulent mixing tests in a rectangular duct and made 
probe measurements of temperature, pressure, and composition within the 
test section. Such measurements are useful to test the present program, 
especially the simplification of equilibrium chemical reactions. 

The test section of the experiment is sketched in figure 19. A high 
pressure gas generator supplied either Mach 2.44 vitiated air at static 
temperature of 1270 K or Mach 2.44 inert gas of 1150 K and a static pressure 
of 0.1 MN/m^, Hydrogen also at 0,1 MN/m^ was injected at Mach 1 (at z = 0) into 
the flow in the heat sink combustion duct from a stepped-wall injector. 

The cross-section of the duct expanded linearly from 5.10 by 9.38 centi- 
meters at z = 0 to 5.10 by 10.48 centimeters at the exit (z = 35.6 cm) to com- 
pensate for boundary layer buildup. The total temperature and composi- 
tion profiles were recorded at the injection step and at the exit 
station; in addition, pitot pressure profiles also recorded at the exit 
station. Detailed test conditions, probe measurements and some 
deduced flow variables like velocity and temperature were reported in 
references 23 and 24. 


26 



In the application of the SHIP program, the velocity and tempera- 
ture from references 23 and 24 at the injection step (2 = 0) were used as 
the initial data which were assumed uniform in the transverse (x) direction. 
Pressure and composition were taken to be uniform at the initial station. 

At the walls, constant wall temperature of 298 K was assumed and heat 
transfer was allowed across the wall boundaries.. Both t^^^o- and three- 
dimensional computations were performed. In the two-dimensional compu- 
tation, the top (north) and bottom (south) boundaries were walls and the 
east and west boundaries were symmetric surfaces. In the three-dimen- 
sional computation, all boundaries were kept the same as the two-dimensional 
computation except that the west boundary (x = 0) was changed to a wall. 

Initial turbulence kinetic energy and dissipation energy rate were 
estimated based on the mixing-length hypothesis. Due to the two-dimensional 
nature of the initial data, two-dimensional turbulent boundary layers were 
assumed adjacent to the walls; the local mixing length was then related to 
the corresponding boundary layer thickness. (See eqs, (4a) and (4b).) in 
the wall jet, two-dimensional velocity and temperature variations were 
assumed and the thickness of boundary layer was taken to be the half of 
the wall-jet height (0.2cm). In the test section, the boundary layer thick- 
ness was assumed to be 0.4 cm. It might be noted that the measurement 
of total temperature profile (in the y-direction) at the injection step 
indicated the variation in a thick layer of 2 cm. It was not clear 
whether the variation in this thick layer was- entirely due to the turbulent 
boundary layer or due to other effects such as wave interactions in the duct. 
Furthermore, no information about flow variations were available in the- 
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transverse direction (the small dimension of the duct) which could help 
to estimate the boundary layer thickness. The estimate of initial mixing 
length and turbulence kinetic energy based on this thick layer (2.0 cm), 
however, was found to result in too much turbulent mixing at downstream 
in both combustion and mixing cases. 

Computations were performed with the values of the following empirical 
constants associated with the "k-e" two equation turbulence model: 

C-j = 1.44, Cg = 1.92, Cp = 0.09, Pr^^^ = 0.7, Pr^ ^ = 1.0 for 0 = u, v, w, and 
k, 1,3 for 4> = e» and 0,5 for ij) = H and f. Computational results at the exit 
station (z = 35.6 cm) were compared with the experimental data in fig. 20 
for the combustion case and in fig. 21 for the pure mixing case. Compari- 
sons show generally good qualitative agreements. However, some discre- 
pancies do exist in these comparisons, for example the composition pro- 
files in fig. 20a. Such discrepancies may be attributed to the simpli- 
fication of equilibrium chemical reactions, the uncertainties associated 
with the initial mixing length and turbulence kinetic energy, and the two- 
dimensional initial data. It has been shown, for example, that nonequili- 
brium chemistry can shift the concentration profiles (ref. 25). Variation 
of initial mixing length could also alter the shape of computational pro- 
files. Furthermore, the estimate of initial turbulence kinetic energy 
by means of the mixing length hypothesis would be significantly different 
if the initial flow fields were three-dimensional. From the geometry of the 
duct, the turbulent boundary layers at two side-walls may contribute sig- 
nificant mixing effects because of the smaller transverse dimension of 
the duct. The two-dimensional flow field at the initial station, however. 
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has resulted in only moderate differences between the three- and two-dimen- 
sional computations. Because of the absence of chemical reactions, the 
composition profiles of the pure mixing case (fig. 21a) are in better 
agreement than those in the combustion case (fig. 20a). In addition, it 
may be noted that combustion builds up pressure in the combustor duct. 

The pressure propagates upstream through the subsonic boundary layer adja- 
cent to the wall and affects the supersonic upstream flow field. Such a 
supersonic-subsonic interaction is not considered in the present para- 
bolic flow program; its effects are difficult to estimate in the present 
comparison. 


DISCUSSION OF THE PROGRAM CAPABILITIY 

The theoretical formulation and numerical procedure of the parabolic 
flow computer program SHIP have been reviewed; In order to achieve econom- 
ical operation of the program, many mathematical and physical simplifica- 
tions were introduced in its development. Some simplifications affect 
only the accuracy of the numerical result; others limit the appli- 
cation of the program. A numerical evaluation of the program has been 
performed for several two- and three-dimensional turbulent, reacting and 
nonreacting flow fields. Generally good numerical predictions are obtained 
when the simplifications on which the program is based are justified. 

The results of application of the present program to three-dimensional 
reacting and nonreacting flow fields are not conclusive because of the 
lack of detailed three-dimensional measurements. However, the program is 
capable of predicting three-dimensional flow fields with qualitatively 
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good agreement with available experimental measurements. Quantitatively 
accurate predictions are, in general, dependent on the proper estimate of 
the initial turbulence kinetic energy and dissipation energy rate and, in 
reacting flows, on the "correct" turbulence-chemical reaction model. For 
chemical reactions in turbulent flow, studies are still needed to define 
a "correct" turbulence-chemical reaction model, and such studies should be 
conducted in a flow field of simpler geometry, e.g. two-dimensional planar 
or axisymmetric flow (ref. 26). With confidence in the turbulence-chemical 
reactions, such a model could be incorporated into the present program. 

The difficulties associated with the estimate of initial turbulence 
kinetic energy and dissipation energy rate are almost proportional to the 
complication of the flow field at the initial station. The complication 
of the initial flow field is usually due to the occurrence of the recircu- 
lating or "elliptic" flow field in the very near field; for example, the 
recirculation immediately downstream of the normal injection, and the 
"elliptic" flow near the trailing edge of a splitter plate. Therefore, 
to reduce the uncertainties of estimating initial turbulence kinetic energy 
and dissipation energy rate, and to make the computational approach 
applicable in the entire flow field, a method for predicting recirculating 
flow field is needed. For the application to supersonic combustor develop- 
ment, the interaction between the supersonic and subsonic flow fields 
becomes important. The downstream pressure effects which are neglected 
in the present parabolic flow computer program should also be considered. 
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CONCLUDIt'IG REMARKS 


The three-dimensional parabolic flow computer program SHIP has 
been evaluated both analytically and numerically. To achieve the effici 
ency in both computer storage and computing time, many mathematical and 
physical simplifications have been introduced into the program. A 
numerical evaluation has been performed for several two- and three-dimen 
sional turbulent, reacting and nonreacting flow fields. Good predic- 
tions are generally obtained from the program when the simplifications 
are justified. 

For the application to supersonic combustor development, the 
present computer program is mainly limited mathematically by the para- 
bolic flow assumption and physically by the equilibrium chemistry simpli 
fi cation. Continuing studies to remove these limitations are needed. 
Provided that recirculation (if present) and downstream pressure effects 
can be properly modeled or separately calculated and the reaction model 
can be improved, the present parabolic flow computer program is capable 
of predicting complicated flow fields in combustors. Thus, the present 
program is potentially valuable for supersonic combustor development and 
experimentation. 



•APPENDIX 

DERIVATION OF THE DIFFERENCE EQUATION 


The derivation of the general form of the finite-difference equa- 
tion (2) from the differential equation (1) is outlined in this Appendix 
for completeness. Details can be found in many publications of Spalding's 
group (for example, refs. 1, 4, 9). 

Equation (1) can be transformed into a finite-difference equation 
by integrating it over the control volume surrounding an arbitrary mode 
P shov/n in figure 22 by dotted lines. W, E, S, N are, respectively, the 
neighboring nodes of P in the west, east, south, and north directions. 

The points, w, e, s, n are the midpoints of the lines PW, PE, PS and PN, 
respectively. In the z-direction , c|) varies in a stepwise manner; this 

makes the finite-difference scheme fully implicit. For the calculation 
of the z-direction convection and of source terms, the variation of 4) 


in the xy plane is also taken to be stepwise. That Is, in the xy plane, 
the value of 4> is assumed to remain uniform and equal to 4> over the 

* r 

dotted rectangle (fig. 22) surrounding the point P and to change sharply 
to 4)^^, 4>^, 4*5 or 4>j^ outside the rectangle. For the cross-stream convec- 
tion from the xz and yz faces of the control volume, the value of 4> 
convected is taken to be the arithmetic mean of the 4> values on either 
side of that face. For diffusion across the xz and yz faces of the control 
volume, 4> varies linearly between two neighboring grid points. 

When the integrations of various terms in equation (1) are expressed 
in the manner described above, the following form of difference equations 
is obtained. 
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F^j = {Ax)(Ay)(pw)p^y/Az 

= (Ay)(pu)^/2 
= (Ax)(pv)y/2 


fp = K H 


(Ay)/5x 


= P. (Ax)/6y 
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Rearranging the terms in eq. (Al) yields equation (;2),, 

♦p ' ¥n ''s'fs ¥e Vw ^ ® («) 

where Aji - (tJI - Lj|)/Ap 
As - (T^ + Lp/Ap 

Ag = (T^ - L^)/Ap (A3) 

\ " '■*)/Ap 

® * <Vp,u ^ 

\ = (tJ; - ■-?;) * <1 * Lp H- (T^ - Lp + + Lp 

+ Fu - Sp (A4) 

When the lateral convection {denoted by the symbol L) is large, some 
of the coefficients Aj^, A^, A^ and Aj^ can become negative, a modifi- 
cation by the "hybrid" scheme is used (ref. 17). The modification con- 
sists of replacing all the T's by T's defined by 

T = i{T + |L| + |T -“iLl|} (A5) 
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Is always positive. 
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F-ig,ure 3^-- A typical comparison of velocity- prafiles at z = 0.487 
based on two different shear stresses in the wall function. 
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Figure 6.- Initial velocity and density profiles of a N 2 -He mixing layer. 
(W^/W^ = p^/p 2 = 7, = -0.0254 m, z =,0.0254 m) 
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Figure 7.- Velocity and density profiles of a N 2 '-He mixing layer. 
(W^/W^ - p.|/p 2 = 7, Zp = -0.0254 m) 
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Figure 8.- Initial velocity and density profiles of a N 2 “He mixing layer. 
(W^/W 2 = 7, p^/P 2 = 1/7, Zj, = 1.27 X 10'^ m, z - 0.0254 m) 
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Figure 9.- Velocity and density profiles of a N 2 "He mixing layer. 
(W^/W 2 = 7, p^/pg = 1/7, Zq = 1.27 X lO"^ m) 
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(b) z = 0.0635 m. 
Figure 9.- Continued 
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Figure 10.- Typical velocity and density profiles of a N 2 -He mixing layer 

with a density-dependent mixing length. 

{W^/W^ = 7, p^/p 2 = 1/7, Zq = 1.27 X 10“^ m, z = 0.0889 m) 
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Figure 13.- Comparison of pressure profiles across the channel. 



Figure 14.- Flow field and injector arrangement of a three-dimensional normal 

injection problem. 
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Figure 15.- Initial profiles and experimental data, 
{s/d = 6.25, (z-z.)/d = 30) 
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Figure 15.- Continued.' 
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Figure 15,- Concluded. 




(a) Velocity. 

Figure 16.- Comparison of computation with experiment, 
(s/d = 6.25, (z-Zj.)/d = 60) 
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Figure 17.- Comparison of computation with experiment, 
(s/d = 6.25, (z-Zj)/d = 120) 
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Figure 18.- Comparison of computation with experiment, 
(s/d = 12.5, {z-Zj)/d = 60) 





y (cm) 



T (°K) 


(b) Temperature. 

( . 

Fiaure 18.- Continued. 



Computation 

iL. 0 

d 

^ = 5.0 

d 

Experiment (ref. 22) 

O ^ = 0 

d 






Figure 19.- Sketch of the test section of a combustor duct. 
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(c) Total temperature ratio profiles = 2380 K). 

Figure 20. >• Continued. 




(d) Pitot pressure, ratio profiles (Pj.gf = 17.1 x 10® N/m^). 


Figure 20.- Concluded. 



(a) Composition profiles. 

Figure 2.1.- Comparison o-f computation with experiment.. Turbulent mixing' of 

hydrogen and inert gas (z = 35.6 cm).." 
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(b) Mach number profiles. 
Figure 21 .- Continued. . 
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(c) Total temperature ratio profiles (T^gj = 2276 K) 
Figure 21.- Continued. 
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(d) Pitot pressure 



Figure 22. The control volume (ref. 4). 
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